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We develop from first principles the coupled wave equations that describe 
polarization-sensitive parametric amplification based on four-wave mixing 
in standard (randomly birefringent) optical fibers. We show that in the 
small-signal case these equations can be solved analytically, and permit us 
to predict the gain experienced by the signal beam as well as its state of 
polarization (SOP) at the fiber output. We find that, independently of its 
initial value, the output SOP of a signal within the parametric gain bandwidth 
is solely determined by the pump SOP. We call this effect of pulling the 
polarization of the signal towards a reference SOP as polarization attraction, 
and such parametric amplifier as the FWM-polarizer. Our theory is valid 
beyond the zero polarization mode dispersion (PMD) limit, and it takes into 
account moderate deviations of the PMD from zero. In particular, our theory 
is capable of analytically predicting the rate of degradation of the efficiency 
of the parametric amplifier which is caused by the detrimental PMD effect. 
©2013 Optical Society of America 

OCIS codes: 230.5440 Polarization-selective devices; 060.4370 Nonlinear 
optics, fibers; 230.1150 All-optical devices; 230.4320 Nonlinear optical devices 

1. Introduction 

Recent years have witnessed a substantial growth of interest in developing nonlinear- 
optical techniques for the control of the state of polarization (SOP) of light beams. The 
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motivation behind such research activities is twofold. First of all, nonlinear optical 
techniques may permit replacing the inefficient and lossy method of polarizing a 
light beam by conventional passive linear polarizers with a the lossless polarization 
attraction of an arbitrary initial SOP towards the desired SOP at the output of a 
nonlinear medium. A key advantage of using lossless polarization attraction is that, 
in contrast with passive linear polarizers, input signal SOP changes do not lead to 
output signal intensity fluctuations or relative intensity noise (RIN). The second goal 
is to flnd efficient ways to excercise all-optical control over the SOP of a signal beam 
by exploiting its nonlinear interaction with a pump beam with a well-determined SOP. 
Here we analyse a novel method for achieving the all-optical control of the SOP of a 
signal beam, namely exploiting the four-wave-mixing-mediated process of parametric 
amplification in a standard telecom optical fiber. 

In short, nonlinear-optical methods allow for designing novel types of polarizers 
with much greater functionality than conventional passive linear polarizers. So far, 
two distinctly different types of nonlinear-optical polarizers were proposed. The first 
class comprises the so-called nonlinear lossless polarizers (NLPs), which are based on 
the cross-polarization modulation (XPolM) of two intense beams in a Kerr medium. 
To the second class belong the so-called Raman polarizers, which are based on 
the polarization-sensitive Raman amplification of a signal beam in a Raman-active 
medium. These two types of polarizers exploit the two complementary manifesta- 
tions of the cubic nonlinearity of fibers - conservative for inducing XPolM effect, 
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and dissipative which is responsible for the Raman effect. Here, we exploit the same 
cubic nonlinearity, more precisely its conservative part, for initiating the process of 
polarization-sensitive four-wave mixing (FWM) between three beams. 

The first NLP was proposed and experimentally demonstrated by Heebner et al 
in Ref. [1]. It was based not on the Kerr nonlinearity, but on a photorefractive ef- 
fect. This polarizer was capable of transforming, in a lossless manner, a light beam 
with an arbitrary initial SOP into a beam with one and the same SOP towards its 
output. The principle of operation of this device was the conversion of energy from 
one polarization component of the beam into its orthogonal polarization component. 
Photorefractive materials are characterized by a nonlinear response which is far too 
slow to be useful in contemporary ultrafast optics. In contrast, the Kerr nonlinearity 
of silica is virtually instantaneous, which makes optical fibers a promising medium for 
implementing lossless polarizers within high-bit-rate telecom networks. The progress 
in developing fiber-based NLPs started from impractical isotropic fibers [2-4] and 
evolved towards cheap and reliable telecom fibers [5-7] or specialty fibers such as 
highly-birefringent and spun fibers [8]. The mathematical aspects of the problem were 
studied in Refs. [9-13], and allowed to get further insight into the physics of fiber- 
based NLPs, whose principle of operation is different from that of photorefractive 
lossless polarizers. Instead of the self-interaction of a single beam in a photorefractive 
material, a two-beam cross-interaction (namely, XPolM) is used in the Kerr medium. 
Namely, an auxiliary pump beam with a well-defined SOP is employed, serving as a 
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polarization reference for the signal beam with arbitrary initial SOP. As previously 
outlined, when using lossless polarizers input signal SOP fluctuations do not lead to 
output RIN [5]. 

Another type of nonlinear-optical polarizer is the Raman polarizer. It is different 
from conventional Raman amplifiers by its sensitivity to the SOP of the pump beam. 
The signal which experiences Raman amplification acquires an SOP which is dictated 
by the SOP of the pump. In this way we may exercise an all-optical control over the 
polarization of the signal beam. Note that conventional fiber-optic Raman amplifiers 
operate in the regime where the output SOP of the signal is independent on the pump 
SOP. The first Raman polarizer was demonstrated by Martinelli et. al. in Ref. [14], 
followed by a number of theoretical papers [15-21]. Similar polarization-sensitive am- 
plification was predicted theoretically and confirmed experimentally in Ref. [22] for 
the Brillouin amplification of a signal beam in standard optical fibers. These devices 
can be similarly called Brillouin polarizers. Since they are based on a gain mechanism, 
which is maximum whenever the signal and pump SOPs are aligned and zero when 
they are orthogonal, in general both Raman and Brillouin-based polarizers suffer from 
severe output RIN in the presence of input signal SOP fluctuations. 

A common feature uniting all of these nonlinear fiber-optic polarizers is that they 
can operate efficiently only in the limit of vanishing polarization mode dispersion 
(PMD). PMD is the effect which is caused by random variations of the magnitude 
and/or orientation of the birefringence along the fiber length, and it is acquired as a 
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result of inevitable technical imperfectnesses in the process of drawing a fiber from 
a preform. Recent progress in fiber manufacturing brought to the market fibers with 
much lower values of PMD than it was previously available. It is this technologi- 
cal breakthrough that made it possible the observation of the previously discussed 
polarization-sensitive effects in optical fibers. Theoretical estimates show that the 
smaller the PMD coefficient, the shorter the total fiber length and the smaller the 
frequency separation of the signal and the pump beams, the better the performance 
of all of the above described polarizers. It is one of the main goals of a theory to be 
able to predict the degradation rate of useful polarization attraction effects which is 
caused by PMD. Such degradation rates for NLPs and Raman polarizers have been 
calculated analytically in Refs. [21,23]. 

It is important to note that the concept of all of these smart polarizers is not limited 
to fiber-optics applications only. Indeed, nonlinear polarizers can be implemented with 
any optical waveguide exhibiting Kerr and/or Raman nonlinearity. Using integrated 
optics waveguides may lift the problems which are associated with fiber PMD, and 
even make nonlinear polarizers less bulky and more compact, providing that the 
waveguide material exhibits nonlinear coefficients which are much larger than silica. 
For example, the silicon-based Raman polarizer proposed in Ref. [24] is free of the 
PMD-induced degradation and has a cm-long size as compared to the km-long fibers, 
thanks to three-four orders of magnitude Raman gain enhancement in silicon with 
respect to silica. 
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The present theoretical study extends the concept of nonhnear polarizers to the 
FWM process in telecom fibers. The goal here is to find the conditions upon which 
the process of parametric amplification is sensitive to the SOP of the pump beam. In 
this way we arrive to the notion of a FWM-polarizer, meaning that the SOP of the 
amplified signal beam is determined by the SOP of the fully polarized pump beam. 
We derive here the coupled wave equations for the pump, idler, and signal beams. In 
the limit of zero PMD, these equations reduce to the equations which were previously 
derived by McKinstrie et. al. in Ref. [25] for describing degenerate FWM in standard 
fibers. The major advantage of our theory is its applicability (slightly) beyond the 
zero-PMD limit, in the sense that it is capable of predicting the degradation rate of the 
efficiency of the FWM-polarizer for low-PMD fibers as well. Knowing this degradation 
rate allows one to properly design practical fiber based nonlinear polarizers. The 
present work substantially extends to the case of random birefringence telecom fibers 
a previous study of polarization attraction in deterministic, high-birefringence optical 
fibers [29]. Note that the polarization-sensitive parametric amplification in optical 
fibers was studied theoretically by Lin and Agrawal [26,27], and also theoretically and 
experimentally by Freitas et. al. in Ref. [28] and resulted in a proposal of a fiber-based 
polarization switch. As discussed in Ref. [29] , FWM-based polarizers are based on the 
polarization sensitivity of parametric gain. Such gain is maximum for a signal SOP 
which is aligned with that of the pump, and zero for a signal SOP orthogonal to the 
pump. Thus FWM-polarizers are not immune from output RIN resulting from input 
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signal polarization fluctuations. Nevertheless, since parametric gain is generally larger 
than Raman gain in silica flbers, FWM-polarizers may employ shorter flbers or lower 
pump powers than Raman polarizers. In addition, since the repolarization capability 
of FWM polarizers is based on parametric gain, these devices provide a more flexible 
control over the gain and repolarization bandwidth. Indeed, such bandwidth may be 
extended up to 70 nm and even include the normal dispersion regime by properly 
engineering the wavelength dependence of the flber dispersion and by adjusting the 
pump power [30]. 

2. Equations of the model 

We shall consider the process of degenerate four- wave-mixing (FWM). This process in- 
volves three waves with frequencies that satisfy the matching condition 2up = Ug + uJi. 
Pump, signal, and idler waves are labeled correspondingly as p, s, i. All three waves 
are co-propagating along the z direction in a telecommunication (i.e. randomly bire- 
fringent) flber. The vectorial theory of parametric ampliflcation in flbers was devel- 
oped in Refs. [25-27], basing on the tensorial properties of silica in the telecom band. 
The starting equation is derived under standard for nonlinear optics approximations, 
from Maxwell's equation with a polarization which takes into account the nonlinear 
cubic response of silica and the birefringence of the flber. Utilizing the Jones rep- 
resentation and transforming into a reference frame moving with the speed of light 
at the corresponding frequency Uf (with / = {p,s,i})), the equations for the Jones 
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vectors of the pump and the signal read as 
AB{ujp, z)Up 



dz 



+\i[{u,-u;)Up + \{Up-Up)u;] 

+\i[{Us ■ u:)Up + (u, ■ Up)u: + (Up ■ u:)Us + (u, ■ u*)Up + (u, ■ Up)u* + (Up ■ u*n] 
+'^jexp{tAkz)[{u, ■ Us)u; + {u, ■ u;)Us + {Us ■ u;)u,] = o , (i) 



s 

s 



\i[{Up ■ u;)Us + {Up ■ Us)u; + (f/, ■ u;)Up + {u, ■ u*)Us + {u, ■ Us)u* + (f/, ■ u*)u,] 

-l^exp{-tAkz)[^{Up ■ Up)U* + {Up ■ U*)Up] = , (2) 



while the idler equation is obtained from Eq. (2) by exchanging labels s and i. 

The Jones vectors Uf = {uxf.UyfY are two-component vectors with Uxf{z,t) and 
Uyf{z, t) being the amplitudes of the polarization components in a fixed laboratory ref- 
erence frame (x, y). Note that the last lines in Eqs. (1) and (2) represent the so-called 
energy-exchange terms. They are responsible for the transfer of energy between differ- 
ent waves, and as such are most important for our analysis of parametric amplification. 
The wave vector mismatch AA; = /3ei+/3es— 2/3ep = Poi+(3os — '^Pop should be as small as 
possible for efficient energy conversion to take place. Here /3e/ and Pof are the propa- 
gation constants of the modes aligned with extraordinary (e) and ordinary (o) axes at 
frequency Uf, respectively. The 2x2 matrix AB{ujf) = A f3{ijjf) {cos 9a3 + sin 9ai) rep- 



resents the birefringence tensor, where A/3(cjj) is the value of the fiber birefringence 
at frequency Uf. Moreover, 6 is the angle of orientation of the axis of the birefringence 
with respect to the fixed reference frame which is defined by the polarization modes 
Bx and By. (T3 and cri are the usual Pauli matrices. 

The orientation angle 6 is randomly varying in fibers used for telecommunication 
applications, which explains the term randomly birefringent fibers that is applied to 
them. In principle, the magnitude of the birefringence AP{ujf) also varies stochasti- 
cally along z. However, as noticed in Ref. [31], the two approaches, one in which 6 
is the only stochastic variable, and the second, where both 9 and A(3 are stochastic 
variables, produce nearly identical results. Thus, here we shall develop our theory by 
assuming the single stochastic variable 6. The angle 6 is driven by a white noise pro- 
cess dzO = ge{z), where {ge{z)) = and {gg{z)gg{z')) = 2L'^^6{z — z'). Here is the 
correlation length which characterizes the typical distance over which 6 changes ran- 
domly. The theory developed below is the natural extension of the one-beam theory 
of Wai and Menyuk in Ref. [31], and two-beam theory of Kozlov, Nuiio and Wabnitz 
in Ref. [7] to the case of three interacting beams. All details of the derivations of the 
final equations of motion with deterministic coefficients starting from Eqs. (1), (2) 
with stochastic coefficients, as well as the approximations which appeared on the way. 
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be found in Appendices A and B. Here we write down the final result: 



' dz + dt 



9 

+ ^ ( 2Ca{z) + =-Cl,{z) ) (PuMU'^"' + ? f2 + ICaiz) 



Xi 



2 

X2 



+ - ( 2 + -Ca(2;) ) 01p01.0ip + + q 2 - 1 VlsV2iV2i 



2 / 1 



+ 3(2- -Ca{z) \ 0ls02p02p + g<^a(2)01p02s02p + g01s02s0L = 
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dz " dt 



X3 



2 / 2 



+ ^ ( 2 + ) 02p02.0;p + ^0L0L + 9(2- JC'al^;) I V2s(FliVli 



2 / 2 \ 3 3 

+ - ( 2 - gCal^) j <p2s<Plp<P\p + -Ca(2:)02p01.01p + -02.01S0L = (5) 

Here Ca{z) = exp[-(3/3)(A(-))2L,^] and Cb{z) = exp[-{l/3){A^-^fLcz], where 
A*^") = A(3{ujp) — A(3{uji). The equation for (pu is obtained from Eq.(4) by exchanging 
the labels i and s; the equations for 02p and 02j are obtained from the equations for 
(pip and (pii, respectively, by exchanging the labels 1 and 2. The polarization compo- 
nents (pif and 02/ are obtained from the original components u^f and Uyf by means 
of a unitary tranformation of the reference frame, see Appendices A and B for details. 

When the value of A*^~) is close to zero, the ^-dependent coefficients Caiz) and Cb{z) 
are both equal to unity, and we restore the model equations derived in Ref. [25] start- 
ing from the Manakov equation. This limit corresponds to vanishing PMD, and it is 
quite natural to call it as Manakov limit. In this limit the conversion of the pump en- 
ergy into the signal is maximally efficient, and therefore parametric amplifiers should 
be designed in such a way that the FWM diffusion length Lf^^ = [(3/3)(A(^))2L,]^i 
is much longer than the fiber length L. To the best of our knowledge, our theory 
for the first time analytically predicts the length scale {Lf^rn) oi degradation of the 
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process of parameteric amplification in telecom fibers. Strictly speaking, our theory 
is valid in two limits: L ^ Lfwm and L ^ Lfwm- In the opposite limit (the limit 
of large PMD, which we call diffusion limit) where L ^ Lf^mi the FWM process 
is totally suppressed. Therefore this regime is not interesting from the viewpoint of 
frequency conversion. Most likely, the intermediate case of L ~ Lf^rn can be ade- 
quately treated only numerically, however we believe that the exponential decay of 
the nonlinear coefficients provides a qualitatively correct description of the rate of 
degradation. 

Another length scale (L^ = [{l/?>){^S-^fL^]-^) appears in Eqs.(3-5), which is called 
the PMD diffusion length. This length is a characteristic which also enters the the- 
ory of two-beam nonlinear interactions. This length scale was introduced by Lin and 
Agrawal in Ref. [23] in the context of fiber-optic Raman amplifiers, and it was identi- 
fied as the typical length at which the mutual orientation of the states of polarization 
of the pump and signal beams is scrambled as a result of PMD. It is quite remarkable 
that this very same length scale not only characterizes the "polarization memory" 
of Raman interactions, but it also characterizes the degradation of cross-polarization 
modulation (XPolM) mediated Kerr-interactions of the two beams, as it was shown 
in Ref. [21]. 

The comparison of the two length scales, Lf^m and L^, shows that polarization 
sensitive frequency conversion is a more demanding process than polarization sen- 
sitive Raman amplification and XPolM-induced polarization attraction in telecom 
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fibers, in the sense that its degradation is characterized by an eight times faster 
spatial degradation rate. In order to overcome this detrimental degradation process 
an experimentalist needs to select an ultra-low PMD highly nonlinear fiber. As the 
degradation rate depends quadratically on the frequency difference between signal 
and pump, the effect of PMD can be also viewed as setting an upper limit to the 
bandwidth of the simultaneous parametric amplification and repolarization process. 

3. Polarization attraction: analytical results 

In this section we shall apply Eqs.(3-5) to describe the effect of polarization attraction 
of a signal or idler wave towards the SOP of a copropagating pump beam by means 
of FWM in a randomly birefringent telecom fiber. We will limit our analysis here to 
the small signal case, i.e., we make the undepleted pump approximation. As we shall 
see, this approximation permits us to obtain relatively simple analytical results for 
the effective bandwidth and gain of the polarization attraction process. From Eq.(3) 
one obtains the two polarization components of the pump amplitude as 

(pip = VPexp{i9ipo + i-f{8/9)Ptotz) 

(f>2p = ^/Q exTp{i92pQ + i-f{8/9)Ptotz) (6) 

where P = |0ip(O)p, Q = |02p(O)P, V^exp(^ipo) and \/Qexp(^2po) are the input 
pump amplitudes in the fiber and Ptot = P + Q is the conserved total pump power. 
Labels X(i 2,3,4) in Eqs.(3-5) indicate the four different FWM processes leading to 
sideband gain through the conversion of two pump photons in two sidebands photons. 
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For vanishing PMD (Manakov limit) A^") = 0, Ca{z) = Cb{z) = 1: peak sideband 
gain is obtained at a frequency detuning Aupi such that the dispersive mismatch is 
compensated by the pump-induced nonhnear phase shift, i.e., (32Au p^ + 'j {16/ 9) Ptot = 
0. In the absence of higher-order dispersion, this condition can only be reached in the 
anomalous dispersion regime (i.e., with < 0). In the opposite case of large PMD 
(diffusion limit) one has Ca{z) = Cb{z) = 0, so that the FWM terms are effectively 
suppressed. Yet, approaching the diffusion limit is equivalent to reducing the effective 
pump power to zero, which correspondingly leads to a peak gain for Aujp2 = 0. Let 
us consider now the intermediate case of L ~ Lfu,m, where peak gain is observed at 
an intermediate sideband detuning Aup^ where Aup2 < Aup-^ < Aupi. In order to 
quantify the sideband gain and evaluate their SOP relative to the pump we need to 
solve Eqs.(3-5). Let us apply the change of variables 0i,2(i,s) = 0i,2(i,s) exp(— u;z/2 + 
iOp^q), where v = (32Au'^ — {16/9)'yPtot- By linearizing Eqs.(4,5) for the sidebands, 
one obtains 



-- = i-M(z) 
dz 9 ^ ' 



where = 0^^, 02i, ^L] > and M(z 



Fa{z)P+Fb{z)Q+v/2 Fc{z)P Fo{z)./PQ Fc{z)^/PQ 

-Fb{z)P -Fa{z)P-Fb{z)Q-v/2 -Fc{z)^ -Fn{z),/PQ 

Fo{z)^/PQ Fc(z)y/PQ Fb{z)P+Fa{z)Q+v/2 Fc{z)Q 

-Fc{z)VPQ -Fd{z)^/PQ -Fc{z)Q -Fb(z)P-Fa{z)Q-v/2 

Fa{z) = 47/3(1 + a(^)/3), Fb{z) = 47/8(1 - C,(z)/3), Fc{z) = 
27/3(Ca(2) + Cb(2)/3); Fd{z) = 8-fCa{z)/9. The solution of Eqs.(7) may be 
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written as 



^{z = L)= exp(fi(L))0(^ = 0) 

where Vl{z) is constructed from M{z) as a Magnus series expansion [33]. Whenever 
the 2r-dependent coefficients F(^a,b,c,d){.z) are slowly- varying over L (i.e., Lfu,m > L) 
we may truncate the expansion after the ffist term 

n{L)=ni{L)= [ M{z)dz = M (7) 

SO that we simply replace -F(a,b,c,d)(^) with their average values 

1 '•^ 



F{A,B,C,D) = Y F{A,B,c,D)iz)dz (8) 
^ Jy=0 

which can be analytically calculated since Ca = {1 — exp{—kaL)) and Cb = k^'^{l — 
exp{-kbL)), where ka = {l/3){A^-^)'^Lc and kb = (8/3)(A(-))^Lc. In the anomalous 
dispersion regime and for sideband frequency detunings Au below a certain cut-off 
value Auc, M has an eigenvalue with positive imaginary part, leading to the effective 
(or average) sideband gain coefficient 

gl = ^^^Pl, (^-A + ^Cl + Cl - SCa + eCaCb 

-\l3lAu:^ - + jPtotAu^ (9) 

From Eq. (9) we obtain the cut-off frequency Auc of the gain band, the peak frequency 
detuning Aup^ and effective gain ge,peak as 
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Aoji = 4c 



27\P2\c^L-^ + 82An^^PtotL 

2 _ 8 / ^'Pl,{m\c'L^' - 2An'^P,,,L) \ 
^''^'''^ "iy 27 ^2C^L-^ + 82 An^-fPtotL J ^ ' 

where c is the speed of hght in vacuum. Let us briefly discuss the role of the differ- 
ent physical parameters in determining the sideband gain and its bandwidth. First 
of all, increasing the birefringence strength An or the fiber length L reduces the 
peak gain coefficient as well as the optimal sideband detuning. In order to study 
the polarization properties of the sidebands we need to consider the eigenvectors of 
M. For any frequency detuning Au within the gain band, let us denote by p the 
eigenvector of M which grows as exp{gez). After a relatively short distance into the 
fiber, we may well approximate the sideband fields as ~ Cpexp{gez), where C is 
the projection or scalar product (which we suppose nonzero for simphcity) of the 
input sidebands polarization vector (f){z = 0) on p. The components of p are such 



that p[l]/p[3] = p[2]/p[4] = a/ P/Q- Idler amplitudes (f)u and 02i correspond to the 
first and third components of 0, respectively. Thus their ratio can be expressed as 



u/(p2i = 'pI^/'M = VP/Q- Since (Pu/hi = (0ii/02^)e*^^-'^«, we obtain that 



tli = lEei9p~^(>Q = "^'P^^ = (11) 

02i VQ <p2p{z = 0) 

A similar treatment can be developed for the signal amplitudes too, which proves the 
polarization attraction of both the signal and the idler to the input polarization of 
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the pump. 

In practice, since for a given sideband frequency detuning the effective gain coefficient 
Qe decreases as the fiber length L grows larger, the corresponding strength of polar- 
ization attraction will be reduced whenever the fiber length approaches Lf^m- As a 
matter of fact, in the diffusion limit Qe = and FWM-induced polarization attraction 
is no longer observed. In the next section we will provide a quantitative description 
of the fiber length dependence of the polarization attraction efficiency. 

4. Polarization attraction: examples 

Let us study the efficiency of polarization attraction as a function of fiber length L, 
hence of PMD. Consider a fiber with the nonlinear coefficient 7 = 11.9 km~^, 
dispersion P2 = —0.5 ps^ km^^, and PMD correlation length Lc = 10 m. The chosen 
parameters are typical for highly nonlinear optical fibers. As well known, the SOP 
of each interacting wave may be represented by means of its corresponding unitary 
dimensionless Stokes vector as Sj = [Sij = S^^ {(l)\j(j)2j + (f>ij(f>2j), = S^^ {i(t)\j(j)2j — 
#1.02,), S^, = 5o"/(l0i.f - I02,f )] (j = i.V.s), where S,, = (|0i,f + |02,f )'^'. The 
input CW pump beam power is set to P^t = S^p = 1 W, and its SOP is defined by 
the Stokes vector 5*^ = [a/cTS, y/OA, a/oT]. We set the input signal power to Ps^in = 
1 mW, whereas the idler is zero at the fiber input, as in typical FWM experiments. We 
compared the numerical solution of Eq. (7) with the analytical solution of Eq. (7). As 
initial condition we employed a set of 10000 input signal SOPs, whose corresponding 
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Stokes vectors are uniformly distributed over the Poincare sphere. Figures 1 and 2 
illustrate the dependence on sideband detuning of the signal gain Qs and its output 
DOP, respectively, for four different values of An (namely, An = 0, 10~^, 1.5 ■ 10~^ 
and An = 1 ■ 10^^), and the fiber length L = 300 m. The signal gain was computed 
as Qs = C^L)"^ \og[Ps^out/ Ps,in\y where Ps,out is the output signal power. The output 
DOP was calculated as discussed in [32] . In Figs. 1 and 2 the curves refer to numerical 
solutions, and the dots to analytical solutions: as can be seen, the first-order term 
of the Magnus expansion provides an excellent approximation of the exact solution. 
Figs. 1 show that, as the birefringence An (or PMD) strength grows larger, the signal 
gain gs is progressively degraded; at the same time, both the peak gain frequency 
detuning Au^ and the cut-off frequency Auc shrink towards zero. In addition Figs. 2 
show that the signal DOP is maximum for sideband frequencies close to peak gain 
values: however the peak DOP rapidly drops from unity as the PMD strength is 
increased (i.e., for An > 10~^). It is interesting to point out that, in contrast with 
the case of the signal, the output DOP of the idler (not shown here) remains close 
to unity throughout the entire gain bandwidth. The increased attraction of the idler 
towards the pump is due to the fact that the idler grows from zero at the fiber input, 
hence its projection on the growing eigenvector p is much larger than for the signal. 

In the second example of Fig. (3) we show the signal DOP as a function of the fiber 
length L, for four different values of the sideband frequency detuning Au = Auj/{2'k) 
(i.e., Au = 0.255 THz, 0.350 THz, 0.365 THz, 0.380 THz): here the PMD value is kept 
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fixed to An = 1.5- 10~^. As it can be seen, for Ai/ = 0.255 THz (which corresponds to 
the peak gain value Aws) the DOP is monotonically increasing with distance and it 
approaches the unit value for L > 500 m. On the other hand, Fig. (3) shows that for 
other values of the sideband detuning the output DOP exhibits a damped oscillating 
behavior and it converges to relatively low values after fiber lengths of the order of 
1 km. 

It is useful to visualize the effectiveness of polarization attraction by means of 
parametric gain or FWM by plotting on the Poincare sphere the end points of the 
Stokes vectors corresponding to either the input or the output distributions of signal 
SOPs, corresponding to the results of Figures 1-3. In Fig. 4 we compare the distribu- 
tion of input signal SOPs, which uniformly covers the sphere ( Fig. 4(a) ), to the 
output signal SOP distribution ( Fig. 4(b) ) from a fiber of length L = 500 m with 
birefringence An = 1.5 ■ 10~^: here the signal detuning is Au = 0.255 THz. These 
parameters correspond to the sideband detuning for peak signal gain (see the dotted 
curve in Fig. 2). As it shown by the solid curve in Fig. 3, the output DOP is as high 
as 0.97, which means a nearly full attraction towards the input pump Stokes vector 
Sp. On the other hand, in Fig. 5 we show the output distributionof signal SOPs when 
the sideband detuning is increased up to Au = 0.350 THz. Fig. 3 shows that the 
output DOP is only 0.73 in this case, which results in a relatively poor polarization 
attraction. It is important to point out that the polarization attraction (to the pump 
SOP) behavior which is illustrated in Figs. 4 and 5 does not depend upon the specific 
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input pump SOP which is selected: indeed, the strength of polarization attraction 
only depends on the pump power level. 

5. Conclusions 

In our study we proposed and analysed a novel type of nonlinear polarizer, exploit- 
ing the degenerate FWM process or parametric optical amplification in a standard 
telecom fiber with randomly varying birefringence. In the FWM-polarizer the SOP 
of the amplified signal (or idler) beam is attracted to the SOP of the copropagating, 
fully polarized pump wave. We have derived the coupled wave equations that describe 
the propagation of the pump, the idler, and the signal in the presence or weak PMD. 
Our model substantially extends previous theory of FWM in optical fibers, since it 
may analytically describe the rate of degradation of FWM efficiency and polarization 
attraction for low-PMD fibers. Knowing the spatial rate of PMD-induced degradation 
permits the proper design of practical nonlinear polarizers based on optical parametric 
amplification in km long nonlinear optical fibers. Polarization attraction and control 
by parametric amplification in fibers is potentially applicable to frequency-conversion 
and phase sensitive amplification devices when combined with polarization-sensitive 
optical processing devices (e.g., an heterodyne receiver). In addition, codirectional 
parametric repolarizers based on low-PMD telecom fibers may be used for compen- 
sating ultrafast input signal SOP fluctuations. Although FWM-based polarizers suffer 
from output RIN, however RIN suppression could be obtained when operating the 
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amplifier in the depleted pump regime, as it occurs with Raman polarizers [34]. 
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A. Appendix: Stochastic theory of parametric amplification 

Our goal is to convert the initial equations for the field (1) and (2) with stochas- 
tic coefficients into corresponding equations with deterministic coefficients. In other 
words, we need to find a way to average the initial equations over the ensemble of 
fibers, which represents all possible realizations of the random fiber birefringence with 
a given statistics. Since both initial and final equations are nonlinear, our procedure 
cannot be done exactly and it will require a number of approximations. Thus, the 
final equations will have a limited range of applicability. 

We use the approach first introduced into the fiber optics theory by Wai and 
Menyuk in Ref. [31]. This approach was formulated for a single beam (or pulse), 
and lead to the derivation of the celebrated Manakov equation and its generaliza- 
tion in the form of the Manakov-PMD equation. An extension of this theory for the 
two-beam configuration was undertaken in Refs. [7,16], and led to the formulation 
of the theoretical basis of XPolM-induced polarization attraction effect in telecom 
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fibers and of Raman polarizers. Here, we need to extend tliis tlieory even furtlier by 
fully taking into account the three interacting beams. Given that all these theories 
have very much in common, we shall omit many repetitions and where appropriate 
we simply refer to prior literature for more details. 

We start with the transformation of field vectors from the laboratory (x, y) frame 
into the local reference frame (1, 2), which is defined by the z-dependent orientation 
of the axis of birefringence: '^f = M{z)Uf, where M{z) is the 2x2 rotation matrix 
defined in Eq. (4) of Ref. [7]. Here \E'/ = ("^Ai/, '?/'2/)- AH terms except one in the field 
equations stay immune to this transformation. The only change is the form of the 



birefringence matrix, which now becomes 



AB{ujf) 



\ 



(12) 



y ±lge -A/3(a;^) j 
The next transformation: = Tf{z)'^f, is aimed at the decoupling of the linear 
portions of the field equations. This goal is reached if the transformation matrix 



( \ 

-al{z) al{z) 



(13) 



obeys the following equation 



.dTp 

dz 



ABp ■Tp = 



(14) 



Matrices T, and Tj are defined in a similar way, with 6i 2 and Ci^2 elements used 
instead of 01,2- The unitarity of this transformation is preserved by requiring that 
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|ai(2;)p + 102(2;)^ = + 1^2(2^)^ = |ci(2;)P + |c2(2;)P = 1- Initial conditions for 

the elements of the Tf{z) matrices are to be determined from the requirement that 
$^ = at z = 0. Thus, ai(0) = 6i(0) = ci(0) = 1 and 02(0) = 62(0) = 02(0) = 0. 

The transformation associated with the Tf{z) matrix brings the equations for three 
fields in the form 

i-gj- + l{N,pm + N^pm + AT,^)/ = . (15) 

As expected, in this reference frame the fields are coupled by nonlinearity only through 
three types of cubic terms: SPolM terms Nspm, XPolM terms Nxpmi and energy ex- 
change terms N^x- The number of these nonlinear terms is very large, and we do not 
provide here their detailed structure. Instead, we refer to Eqs. (9)-(12) in Ref. [7] 
where the SPolM and XPolM nonlinear terms are written down explicitely. In our 
present theory we have all these terms as well, and in addition get energy-exchange 
terms in the form of cubic products involving three different fields. 

Coefficients prior to these terms are some self- and cross- fourth-order polynomials 
composed of 01^2(^)5 ^1,2 (2^), Ci^2(^) and their complex conjugates. It is convenient to 
work with quadratic coefficients Um and u*^ (m = 1 -^ 30). Coefficients Um with m = 
1 -i- 14 are identical to those introduced immediately below Eq. (12) in Ref. [7]. They 
are divided into self-terms: ui = |aip — |a2p, U2 = —(0102 + 0*02), = z(oi02 — 0^03), 
U4 = 201O2, = a\ — a*2 ^ Uq = —i{al + 02^), and cross-terms: Uj = 0^61 — O262! 
Us = — (&1O2 + b^al), ug = i(6i02 — 0^63), mio = — ?(o^6i + 02&2)) = 01&2 + bia^, 



29 



Ui2 = aibi — a^bl, M13 = —i{aibi + a2bl), = 2(0162 — a2&i). The three-beam theory 
additionally brings 16 new coefficients. Coefficients Um with m = 15 22 are the 
same as Um with m = 7 ^ 14 but with 61 2 replaced with ci^2- Coefficients Um with 
m = 23 -i-30 are the same as Um with m = 7-i- 14 where 01^2 is replaced with bi^2, and 
simultaneously 61^2 is replaced with ci_2- 

Nonlinear coefficients in Eq. (14) are products of the type UmUn or Umu"^- They 
are 2;-dependent random coefficients, because they depend on the stochastic variable 
ge{z). We need to find average values of all nonlinear terms, which are of the form, 
for instance, ul(j)is(f>2p4>2s- '^^^^ place where the most important approximation 

comes into play. We assume that the following factorization is valid: {ul(f>is4>2p<p2s) ~ 
(ug) (0is02p02s)- This factorization is justified whenever the spatial evolution of the 
fields is much slower than the spatial evolution of the nonlinear coefficients, or vice 
versa, whenever the spatial evolution of the fields is much faster than the spatial 
evolution of the nonlinear coefficients. In the context of parametric amplification and 
in the absence of group- velocity dispersion, the nonlinear evolution of the fields scales 
with the nonlinear length Lj^f^ = ^Pp-, where Pp is the pump power. In its turn, 
the 2;-dependence of (umUn) or {umuD is governed by two different length scales. On 
the one hand we have the relatively short spatial scales which are associated with 
the correlation length Lc and the beat length Lb, both of which are typically less 
than 100 m. On the other hand we have the relatively long spatial scales which are 
associated with the PMD diffusion length and the FWM diffusion length Lf^^. 
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For practically interesting situations we need to provide the following hierarchy of 
scales: Lc, Lb ^ L, L^l <^ La, Lfwm- In this range, the factorization approximation 
is well justified: with this limitation in mind, we may proceed further. 

In order to find averages of the type and |MmP, it is convenient to group 
coefficients as Gi = {ui, U2, M3}, G2 = {u^, M5, Uq}, = {uj, Ug, Ug, Uio}, 

G4 = {Uu, U12, Ui3, Uu}, Gs = {Ui5, M16, Ul7, ^is}, Ge = {Ul9, U20, U21, U22}, 

G7 = {u23, U2A, U25, U26}, and Gs = {u27, U2S, U29, ^3o}- For each group we were 
able to formulate a closed system of linear first-order differential equations by using 
Eq. (14). For an example of such a system, we may refer to Eq. (13) in Ref. [7]. 

Next we need to know the average values of quadratic forms composed by these 
coefficients. They can be found from the solutions to the equations of motion for the 
average of the generic function F. For instance, for F{ui,U2,U3,6) we need to solve 
the equation dz{F) = {G{F)). The generator G is to be constructed by a procedure 
described in the Appendix of Ref. [31]. For a specific example of G, we may refer to 
Eqs. (14) and (20) in Ref. [7]. Note also that the average over different realizations of 
the fiber birefringence can be replaced by a spatial average as 



by assuming that the ergodicity hypothesis is valid. 

With this procedure at hand, we are able to find the mean values of u^, with m = 
9, 10, 13, 14, 17, 18, 21, 22, 25, 26, 29, 30 by solving the equation d,VA = MaVa for 
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(16) 



the vector Va = {{Si), {Si), {Si), {Sj), {S^S,), {S^S,)f, where {S,, S2, S-s, S,} is 
any of the groups Gi with i = 4-^ 8, and with the matrix Ma given by 



c 


c 











2A- 












-2A+ 

















2A(+) 




















-2A( 





A(+) 




























(17) 



Here A(+) = Al3{up) + A(3{uji) and A(~) = A/3(a;p) - A(3{uji). It is a straightforward 
calculation to get an estimate A''~''/A^+) ~ Atu/cOp, where Aco = Up — Ui = Ug — cOp- 
For typical fiber parameters, the evolution associated with A*^"^-* is very fast, while 
A*^") defines a much slower spatial scale. Setting A^~^ to zero defines the Manakov 
limit, and brings us back to the formulation of coupled wave equations with constant 
in z nonlinear coefficients. The difference of A^~^ from zero means the inclusion of 
effects caused by the PMD. In this case we are dealing with 2;-dependent nonlinear 
coefficients. 

Next we calculate the averages of the type [wmP, with m = 
9, 10, 13, 14, 17, 18, 21, 22, 25, 26, 29, 30. Thereto we formulate 

the equation of motion dzVs = MbVb for the vector Vb = 
{{\Si\'), {\S2\'), {\Ss\'), {\S,\'), {S2S*), {S;S,), {S^SD, {SlS,)f where 
{ 5*1, 5*2, 5*3, 5*4} is any of the groups Gi, with i = 4-^8, and where the ma- 
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trix Mb reads as 
/ 



2L- 



2L 



-1 



-2L: 



A(-) A(-) 



-A+ -A- 



A(+) A(+) 







-A(-) -A(^) 



-A(- 



A(+) -A(+) -L-i 



A(+) -A(+) -LZ^ 



A(-) -L-i 



-A(-) A(-) 



-L 



-1 



Note that initial conditions for the averages of the type {um{z)un{z)) and 
{um{z)u*^{z)) can be found from the initial conditions for the coefficients ai^2, 
hi^2i and 01^2, and by observing that (Mm(0)M„(0)) = Mm(0)M„(0). Thus we find 
"1,5,7,12,15,20,25,28(0) = 1 and ^6,10,13,17,21,26,29(0) = while the remaining coefficients 



are all zero. 



Next, we turn to cross-terms like (umUn) with m ^ n. Many of these terms 
are zero, mainly because of the imposed zero initial conditions. Nonzero co- 
efficients are (^14^*22), (^^14^2)5 {uioUis), (niou^g), (^6^29), {ueulg), (^3^25), 
and {usul^). The first four of these coefficients can be found by solving 
the equation dzYs = MbVb with the matrix Mb defined as in Eq. (18), 
and where the vector Vb is identified with any of the following vectors: 

((niiUig), (ni2U2o), (^il3^^2l), {u uU 22) , {U13U20), {U12U21), {uuUig), (^11^22))^, 
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{{unulg), {ui2U*2q), (M13M21), («12^l), {uuulg), {UiiU*22)V, 

{{U7U15), {uguw), {ugUn), (uioUis), {ugUie), {ugUn), (uiqUis), {uyUis))^, 

{{^rul^), (ugule), (ugulj), (uio^g), (^9^*5), (nsn*^), (^10^*5), (uy^g))^. 

In turn, the coefficients (M6M29), {ueuld)^ ('^3'"25), and (M3M25) can be found from 
the equation c^^V^ = MqVc with the matrix Mc defined as 



A(+) A(+) 



-A(+) -L- 



-A(+) 



A(+) 



-L-^ A(- 















2L 



-1 



-A(+) -A(+) 



-2L 



-1 



v 



(19) 



-A(-) -L-^ 

when we associate the vector Vc with any of the follow- 



ing 



vectors: 



{{u6U2g), {U5U29), {ueU28), {u^U2s) , (^4^27), (^^4%o))' 



((M6M29), (M5M29)' (W6M28)' (^5^28)' (W4M27)' (w4M3o))^) (("3M25), (M3M24), (M2M25), (M2M24), (M1M23), 

and ((msM^s), (M3M24)) (W2M25)' (w2M24)> (^i^s), (^i^e))^- 



B. Appendix: Analytic estimation of the nonlinear coefficients 



In this Appendix we look for approximate analytical solutions to the linear systems 
of equations for the vectors Va, Vb and Vc- This task is equivalent to finding the 
eigenvalues and eigenvectors of the matrices Ma,b,c- Additionally, we need to find 
the decomposition of the intial vectors Va,b,c(0) in the basis of the corresponding 
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eigenvectors. In this way we may determine the z-dependence of the nonhnear coef- 
ficients. 

We shall give a detailed analysis for the Ma matrix, and sketch only briefly the 
results for the other matrices. We develop a perturbative approach, by assuming 
that A*^"^ is much smaller than A^"*"^ and First, setting A*^^^ to zero we get 
a much simpler matrix Ma- The difference AM4 = Ma — Ma is therefore a small 
correction. The matrix Ma has a doubly degenerate eigenvalue A^i = and two 
corresponding eigenvectors e^i = (0,0,0,1,0,0)"^ and 6^12 = (1,1,1,0,0,0)"^. The 
other eigenvalues of Ma all have relatively large negative real parts, in the sense 
that the corresponding eigenvectors vanish with distance very quickly. The spatial 
scale of this decay is determined by the correlation length Lc and the beat length 
Lbi both of which are typically less than 100 m. So, the characteristic decay rate is 
estimated as Ltransient ~ 100 m. After the transient decay is over we can write the 
solution of dzVA = MaVa as Va(-z) = (C'iCai + C2eA2) exp{XAz) = CiCai + (^26^2, 
where Ci = (Va(0) ■ e*Ai)/{eAi ■ e^i) and C2 = {Va{0) ■ e*^2)/i^A2 ■ 6^2)' thanks to the 
orthogonality of the set of eigenvectors of Ma- 

When A*^") is different from zero, the degeneracy is lifted and the doubly degenerate 
eigenvalue split into two different eigenvalues A^^i = and Xa2 7^ 0. Let us find 
Aa2 by way of developing the perturbative analysis. First we find the eigenvalue 
equation for the exact Ma matrix. It is det(MA - A/) = 32[A^~^^[A{+)]'^L^^\ + 
(16[A(-)]2[A(+)]2 + 12[A(-)]2L-2 + 12[A(+)]2L-2)A2 + (16[A(-)]2L-i + 16[Ai+)fL-' + 
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4L;3)A3 + (4[A(-)]2 + 4[A(+)]2 + 9L~'^)\^ + 6^"^^ + = 0, where / is the unity 
matrix. Since A^"^ is small, we expect that the correction to the unperturbed zero 
eigenvalue is also small. By keeping in the eigenvalue equation terms no higher 
than second order in A, we get after some simplifications the approximated solution 
Aa2 ~ — (8/3)(A'^~))^Lc = —L'j^^. The eigenvector corresponding to eigenvalue A^n 
(Aa2) is eAi = cai + eA2 (eA2)- The perturbed solution is Va{z) = CiCai exp(AAi2) + 
C2eA2exp(AA22;) = CiCai + (^26^2 exp(-2;/L/^m). Here Ci = {Va{0) ■ e*^^) /{cai ■ e\^). 
The exact expression for 6*26^2 is cumbersome, however under the condition A*^^-* <^ 
A^^\ we can use the equality of Va and Va in the limit of A^"^ — )■ 0, and write 

C2eA2 = CiCAl + 6*26^2 — CiCai- 

Now we can turn to evaluation of the nonlinear coefficients. Let us start with the 
averages (ug) and {uIq). Coefficients Ug and Uiq belong to the group of coefficients 
denoted earlier as G3. For this group, the vector Va contains (wg) and (miq) as the 
third and the fourth element, respectively. The initial condition reads as Va(0) = 
(1,0,0, -1,0, i)^. Thus, for L > Ltransient we find (wg) = (l/3)exp(-2/L/^m) and 
("10) = -exp(-2/Lj^^). Similarly, we find {u\j) = {ul^) = {uD and (uj^) = (w^g) = 
(ulo). With initial conditions Va(0) = (0, 1, —1, 0, —i, 0)^ we get (^^3) = (^^4) = 
{uli) = {uD = {ulg) = {uIq) = 0. 

The matrix Mb can be considered similarly. Again, in the limit A*^"-* — )■ this 
matrix possesses a doubly degenerate zero eigenvalue Xb, while the other eigen- 
values have large negative real parts, so that the corresponding eigenvectors van- 
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ish after a certain propagating distance, say Ltransient- Whenever A*^"^ is differ- 
ent from zero, the degeneracy is hfted and tlie doubly degenerate eigenvalue 
is split into A^i = and Xb2 = —Lfwm- Thus, for L > Ltransient we find 
(l^gp) = (l^irp) = (Kl) = (l/2)-(l/6)exp(-./L;_); {\u,o?) = = 
(l^^eP) = (l/2) + (l/2)exp(-^/L;_); {\u,,\') = {\u,,\') = {\ul,\) = (1/2) + 
(l/6)exp(-z/Ly^„); (|ni4p) = {\ul^\) = (|n3oP) = (1/2) - (1/2) exp(-z/L/^ J. 
By using the same eigenvectors and eigenvalues of matrix Mb we find also that 
(M14M22) = 0, (M14M22) = (1/2) - (1/2) exp(-2;/L/^m), (miqMis) = - exp{-z/Lf^rn), 
and (uio^s) = (V2) + (1/2) exp{-z/ Lf^rn)- 

Finally, matrix Mc possesses a nondegenerate eigenvalue Ac = in the limit A*^~) = 
0, with the other eigenvalues vanishing for z > Ltransient- The perturbative approach 
yields the correction to the zero eigenvalue: Ac = —L'^^. Then, for L > Ltransient we 
find {UQU29) = 0, {uqu*2q) = (2/3) exp(-2/Lrf), (M3M25) = (M3M25) = (1/3) exp(-2/Ld). 

When all these nonlinear coefficients are substituted in the equations for the field, 
we arrive to the final result: Eqs. (3)- (5), which represent equations with deterministic 
coefficients, as desired. 
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ure Captions 

i Figure 1: Dependence of signal gain on its frequency detuning from the pump, 
with L = 300 m. Curves and circles were obtained with z-varying or average 
M coefficients, respectively. Moreover An = (solid curve); An = 1.0 ■ 10"^ 
(dashed curve); An = 1.5 ■ 10^^ (dotted curve); An = 1.0 ■ 10^^ (dash-dotted 
curve) . 

i Figure 2: Same as Fig.l, but for the signal DOP. 

• Figure 3: Signal DOP versus fiber length L with An = 1.5 ■ 10~^, and different 
values of the sideband detuning frequency: Au = 0.255 THz (solid curve); Au = 
0.350 THz (dashed curve); Au = 0.365 THz (dotted curve); Au = 0.380 THz 
(dash-dotted curve). 

i Figure 4: Tips of input (a) and output (b) signal Stokes vectors on the Poincare 
sphere for a fiber length L = 500 m. An = 1.5 ■ 10~^, and A^ = 0.255 THz. For 
the sake of clarity, only 225 vectors are represented instead of the 10000 used 
in the simulations. Input vectors are distributed uniformly over the sphere. The 
empty triangle represents the input pump Stokes vector. 

i Figure 5: Output signal Stokes vectors with L = 500 m and Au = 0.350 THz. 
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Fig. 1. Dependence of signal gain gs on its frequency detuning from the pump, 
with L = 300 m. Curves and circles were obtained with 2;-varying or average 
M coefficients, respectively. Moreover An = (solid curve); An = 1.0 ■ 10~^ 
(dashed curve); An = 1.5 ■ 10"^ (dotted curve); An = 1.0 ■ 10~^ (dash-dotted 
curve). 
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Fig. 2. Same as Fig.l, but for the signal DOP 




L[m] 

Fig. 3. Signal DOP versus fiber length L with An = 1.5 ■ 10~^, and different 
values of the sideband detuning frequency: Az/ = 0.255 THz (solid curve); 
Au = 0.350 THz (dashed curve); Az/ = 0.365 THz (dotted curve); Au = 
0.380 THz (dash-dotted curve). 
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Fig. 4. Tips of input (a) and output (b) signal Stokes vectors on the Poincare 
sphere for a fiber length L = 500 m, An = 1.5 • 10~^, and = 0.255 THz. For 
the sake of clarity, only 225 vectors are represented instead of the 10000 used 
in the simulations. Input vectors are distributed uniformly over the sphere. 
The empty triangle represents the input pump Stokes vector. 
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Fig. 5. Output signal Stokes vectors with L = 500 m and Au = 0.350 THz. 
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